{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Using a numerical integration within fsolve\n",
    "## In this example, we find the chemical potential of a collection of non-interacting identical Bosons.\n",
    "*May 5, 2021*\n",
    "\n",
    "- First import some required modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import quad\n",
    "from scipy.optimize import fsolve\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The left-hand side of the equation that we will solve is:\n",
    "\n",
    "$2\\pi^2 n\\left(\\dfrac{\\hbar^2}{2mk_\\mathrm{B}T}\\right)^{3/2}$\n",
    "\n",
    "where\n",
    "> $n$ is the particle number density<br>\n",
    "$m$ is the mass of the particles<br>\n",
    "$T$ is temperature\n",
    "\n",
    "- We will enter the appropriate values of $^4$He and use a temperature of 4 K."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7937834051194786"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 2.2e28 # 1/m^3\n",
    "m = 6.68e-27 # kg\n",
    "hbar = 1.05e-34 # Js\n",
    "kB = 1.38e-23 # J/K\n",
    "T = 4 # K\n",
    "LHS = 2*np.pi**2*n*(hbar**2/(2*m*kB*T))**(3/2)\n",
    "LHS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the right-hand side of the equation that we want to solve is the integral:\n",
    "\n",
    "$\\int_0^\\infty \\dfrac{x^2}{e^{x^2-\\eta}-1}\\, dx$\n",
    "\n",
    "where $\\eta=\\mu/k_\\mathrm{B}T$ and $\\mu$ is the chemical potential.\n",
    "\n",
    "\n",
    "- First, define a function that represents the integrand of the integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrand(x, eta):\n",
    "    return x**2/(np.exp(x**2 - eta) - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Next, define a second fucntion that calls the first function, integrates it numerically using *quad()* and then returns LHS - RHS where the RHS is the integral.  Ultimately, we want to find the value of $\\eta$ that makes this difference zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(eta):\n",
    "    y, err = quad(integrand, 0, np.inf, args=(eta,))\n",
    "    return LHS - y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Now, use *fsolve()* to find the solution for $\\eta$ when $T=4$~K.  The second option in *fsolve()* is an initial guess for the value of the solution.  We can then use the $\\eta$ solution to calculate $\\mu/k_\\mathrm{B} = \\eta T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value of eta is -0.06714819569069219 when the temperature is 4 K.\n",
      "The value of mu/kB is -0.26859278276276877 K when the temperature is 4 K.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-11-abff491f9b2d>:2: RuntimeWarning: overflow encountered in exp\n",
      "  return x**2/(np.exp(x**2 - eta) - 1)\n"
     ]
    }
   ],
   "source": [
    "sol = fsolve(func, -0.07)\n",
    "print('The value of eta is', sol[0], 'when the temperature is', T,'K.')\n",
    "print('The value of mu/kB is', sol[0]*T, 'K when the temperature is', T,'K.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can put the *fsolve()* statement inside a for loop so that $\\eta$ and $\\mu/k_\\mathrm{B}$ can be found at many different temperatures.\n",
    "\n",
    "- First, delete the value assigned to $T$ and define another function for the LHS of the equation that we want to solve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7937834051194786"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "del T\n",
    "def LHS(T):\n",
    "    return 2*np.pi**2*n*(hbar**2/(2*m*kB*T))**(3/2)\n",
    "LHS(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We also need to update the function \"func(eta)\" so that \"LHS\" $\\to$ \"LHS(T)\" is given as a function of T."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(eta):\n",
    "    y, err = quad(integrand, 0, np.inf, args=(eta,))\n",
    "    return LHS(T) - y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to provide *fsolve()* with initial guesses for each iteration of the loop.  The strategy will be to use the $\\eta$ solution from the previous iteration as the guess for the current iteration.  We then only need to have a reasonable first guess for the first iteration.\n",
    "\n",
    "We will start the loop at the highest temperature when the chemical potential should have a value that is close to the classical result for an ideal gas. In this case, $\\eta=\\mu/k_\\mathrm{B}T$ is expected to be close to $\\ln\\left(n/n_\\mathrm{Q}\\right)$, where:\n",
    "\n",
    "$n_\\mathrm{Q}=\\left(\\dfrac{mk_\\mathrm{B}T}{2\\pi\\hbar^2}\\right)^{3/2}$\n",
    "\n",
    "- Create some empty lists and then execute the for loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-11-abff491f9b2d>:2: RuntimeWarning: overflow encountered in exp\n",
      "  return x**2/(np.exp(x**2 - eta) - 1)\n"
     ]
    }
   ],
   "source": [
    "TList = []\n",
    "etaList = []\n",
    "muList = []\n",
    "T0 = 100\n",
    "nQ = (m*kB*T0/(2*np.pi*hbar**2))**(3/2)\n",
    "guess = np.log(n/nQ)\n",
    "for T in np.arange(100, 4, -0.1):\n",
    "    TList += [T]\n",
    "    sol = fsolve(func, guess)[0]\n",
    "    etaList += [sol]\n",
    "    muList += [sol*T]\n",
    "    guess = sol # Update the guess value before starting the next iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Here's a plot of $\\eta$ as a function of temperature with the temperature axis on a log-scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ20lEQVR4nO3de7zVc77H8ddn1+5y0JTagy6UE7l137pgSoqQUSYpl6GpNCXMjEuYMmcm1yE6mgq5ZUak0KASZVBSaXeR6OjgNFOZYcsccptqfM8fn3qcVKP2Xpfv+q31fj4e66G12uu3PtF6P76+v+/n+7UQAiIiklxFsQsQEZHUKMhFRBJOQS4iknAKchGRhFOQi4gknIJcRCThqsb40Hr16oXGjRvH+GgRkcRaunTpxyGEkp1fjxLkjRs3pqysLMZHi4gklpn9eXevp2VqxcxONbN3zOxdM7s2HdcUEZG9k3KQm1kVYDxwGnAUcK6ZHZXqdUVEZO+kY0TeDng3hPB+CGEzMAXomYbriojIXkhHkDcA1u3wfP2210REJAvSEeS2m9d22YnLzAabWZmZlZWXl6fhY0VEBNIT5OuBRjs8bwh8sPMPhRAmhhBKQwilJSW7rJ4REZFKSkeQLwEOM7MmZlYN6Ac8k4br7mrFCnjlFdiyJSOXFxFJopSDPISwFbgUeB5YDUwNIbyV6nV3a8wYOPFEKCmBvn3h978HTdOISIGzGAdLlJaWhko1BH32GcyZAzNnwqxZ8OGHYAbt2kGPHtCrFxxzjL8mIpJnzGxpCKF0l9cTFeQ7+uYbWL7cQ33mTFiyBEKAI4/00fo55/ivRUTyRP4F+c4+/BCeegoefxzmzfNQb97cA71fP2jaNL2fJyKSZf8qyPNn98MDDoChQ+Hll2H9ehg7FmrVguuvh8MOgxNOgAcegE2bYlcqIpJW+RPkO6pfHy67DF59Ff7yF7j1Vvj4Yxg0CA48EC68EF56yadnREQSLj+DfEeNGsE118Dq1bBwIVxwATz9NJx0ko/Ub78dNm6MXaWISKXlf5BvZwYdOsC998Lf/gaTJ0PDhjB8ODRoAP37w+uvx65SRKTCCifId1SzJpx3njcXrVwJAwbAk09C+/Zw7LHw0EPw1VexqxQR2SuFGeQ7at4cJkyADRtg3Dj48ksP9gYNYMQI+OtfY1coIvKdFOTb1aoFw4bBqlW+8uXEE+GWW6BxYw/2VasiFygisnsK8p2ZQefOviZ9zRpf6TJlio/cTzsN5s71NeoiIjlCQf5dmjaF8eNh3Tq48UbvJD35ZGjTBp54QssXRSQnKMj3Rt26Pl++di3cf7/Po/fp4/u6PPIIbN0au0IRKWAK8oqoUQMGDoS33/bplipV4Mc/hmbNPOA3b45doYgUIAV5ZVSp4htzvfEG/PGPsP/+cPHFPhUzbpyWLopIVinIU1FUBD17eiPR7Nlw8MG+NcC//7vPrf/jH7ErFJECoCBPBzPo3h3mz/c9XJo2hUsv9SmXBx/UHLqIZJSCPJ3MfP35K6/A88/D97/vc+pHHQWPPqpVLiKSEQryTDCDU06BxYt9g66aNeH886FlS5g+XevQRSStFOSZZAZnnunrz6dM8VUtP/qR7+cyZ07s6kQkTyjIs6GoyFe5vPWWb8j18cc+Yu/e3Ve+iIikQEGeTVWr+na577wDd97p54y2bu2vrVsXuzoRSSgFeQzVq8MvfgHvvQdXXeXTLocfDtdeC59+Grs6EUkYBXlMderAbbf5CL1PH/jtb30N+l13qUtURPaagjwXHHII/P73sGwZtGoFP/85HHkkTJumFS4iskcK8lzSurWvZpk9G/bZB845x7fUXb48dmUiksMU5Llme5fo8uV+vujq1dC2re/l8uGHsasTkRykIM9VVarA4MHw3//tN0YnTYLDDoPRozV/LiLfoiDPdbVrwx13+FFznTrB1VfD0UfDM89o/lxEAAV5cjRrBjNm+Px5cbHvuti9uzcZiUhBU5AnzfZu0Lvu8oaili19lYvWn4sULAV5EhUXw+WXw7vv+k3QsWN9xP7II5puESlACvIkq1sX7r7bD7Y45BA/dq5zZ3jzzdiViUgWKcjzQWkpLFwI993n54m2bg1XXAGffRa7MhHJAgV5vigqgkGDvN1/0CD4z//06ZbJkzXdIpLnFOT5pm5duOcen25p1AguuMBPLVq1KnZlIpIhCvJ8VVoKixbBxIke4q1awZVXwuefx65MRNJMQZ7Piop8VcuaNX526Jgxfn7oM8/ErkxE0khBXgjq1vV9WxYs8E7Rnj39yLn162NXJiJpoCAvJB07wtKlvu/57Nm+Ve7YsfDPf8auTERSkFKQm1kfM3vLzL4xs9J0FSUZVFwMw4d7a/8JJ8DPfgbt23vAi0gipToiXwX8CJiXhlokm5o0gVmz4PHHYcMGaNfOd1nctCl2ZSJSQSkFeQhhdQjhnXQVI1lm5odXrF4NP/2p799y1FHw9NOxKxORCsjaHLmZDTazMjMrKy8vz9bHyt6oXRsmTPCboXXqQK9e/tDNUJFE2GOQm9lcM1u1m0fPinxQCGFiCKE0hFBaUlJS+Yolc3a8GfrCCz46v+ce+Oab2JWJyHfYY5CHELqFEI7ZzUP//52Ptt8MXbUKjj0Whg6FLl18LbqI5CQtP5TdO/RQmDsXHnjA9z9v0cJH6lu3xq5MRHaS6vLDs8xsPdARmGlmz6enLMkJZjBggN8M7dEDrr3WV7csXx67MhHZQaqrVqaHEBqGEKqHEA4IIXRPV2GSQw46CJ58Ep54Aj74wKdcrrsOvvoqdmUigqZWpCJ69/bR+UUXwa23+kZc8+fHrkqk4CnIpWLq1PF58zlzYMsW6NQJLrlEh1iIRKQgl8rp1s2PlLviCt+Q6+ijYebM2FWJFCQFuVTePvvAHXf4MXO1a8MZZ8D558PGjbErEykoCnJJXbt23kj061/DtGneSDR9euyqRAqGglzSo1o1+I//gLIyaNDA9zs/7zyNzkWyQEEu6dWiBSxeDKNG+XJFjc5FMk5BLulXXAzXX//t0fm558LHH8euTCQvKcglc7aPzm+4wRuKjj5ao3ORDFCQS2YVF8PIkT46b9hQo3ORDFCQS3a0aAGLFn17dP7UU7GrEskLCnLJnu2j86VLfXTeu7dG5yJpoCCX7Gve3EfnN96o0blIGijIJY7iYhgx4tuj8/PPh08+iV2ZSOIoyCWu7aPzUaNg6lR/Pnt27KpEEkVBLvFtX3e+eLHvrnjaaTBkCHz+eezKRBJBQS65o00bX6Z49dUwcaKvdNF+5yJ7pCCX3FKjBtx2G8yb50fNde7swf7117ErE8lZCnLJTSec4Ic+//SnMHo0lJbCsmWxqxLJSQpyyV377gt33+03P//+d2jf3m+KbtkSuzKRnKIgl9zXvTusWgV9+/pWuccd52eHigigIJekqFMHHnnED674n/+B1q1hzBj45pvYlYlEpyCXZDn7bB+dn3KKnxd60kmwdm3sqkSiUpBL8hx4IDz9NDz0kN8Abd4c7r8fQohdmUgUCnJJJjPo3x/efBOOPRYuvhh++EP4619jVyaSdQpySbZDDoG5c+Guu+DFF310rsMrpMAoyCX5iorg8sth+XJo3NgPrxgwADZtil2ZSFYoyCV/HHEEvPaa76r48MPQsiUsWBC7KpGMU5BLfqlWzfc5nzfPn3fq5MG+eXPcukQySEEu+en4473Fv39/uPlm6NhRTUSStxTkkr/22w8eeMBPH/rzn313xXHjtExR8o6CXPLfWWf5MsUuXeCyy3y/8w8+iF2VSNooyKUwHHQQzJwJ48f7/Hnz5n5eqEgeUJBL4TCDSy7xZYqHHurt/v37w2efxa5MJCUKcik8zZr5MsWRI+EPf/BlijqJSBJMQS6FqbgYbrjBA7yoyE8iuu46LVOURFKQS2E77jhYscI7QW+9FTp0gLffjl2VSIUoyEX22893T5w+Hdatg7Zt4Xe/017nkhgpBbmZ3W5m/2VmK81supnVTlNdItnXq9f/L1O8/HItU5TESHVEPgc4JoTQAlgDXJd6SSIRHXigL1OcMMHnz1u08L3PRXJYSkEeQnghhLB129NFQMPUSxKJzAyGDoWlS6FRIx+pDxkCX34ZuzKR3UrnHPkA4Ll/9ZtmNtjMysysrLy8PI0fK5IhRx4JixbB1VfDvfd6i//y5bGrEtnFHoPczOaa2ardPHru8DMjgK3A5H91nRDCxBBCaQihtKSkJD3Vi2Ra9epw221+eMWmTdC+PYwerRuhklOq7ukHQgjdvuv3zewi4AygawjajUjyVNeusHKlHyl39dUwe7bved6gQezKRFJetXIqcA1wZghBE4iS3+rW9f1Z7rsPFi70G6E6Vk5yQKpz5OOA/YA5ZrbCzO5JQ00iucsMBg2CZcugSRM/Vm7wYPjii9iVSQFLddVK0xBCoxBCq22PIekqTCSnbd+v5ZprvJmoTRtf5SISgTo7RSqrWjVv63/xRR+Rd+zoN0Z1I1SyTEEukqouXfxG6Jln+gi9WzdYvz52VVJAFOQi6bD//jBtmh8t9/rrfiNUB1dIlijIRdLFzHdRXL4cmjb1gysGDYLPP49dmeQ5BblIuh12GCxYAL/8JTz4oN8IXbIkdlWSxxTkIplQXAw33QQvvQRff+37nt9yC/zzn7ErkzykIBfJpM6d4Y03fL35L3/pHaLr1sWuSvKMglwk0+rUgSlTYNIkX2veogVMnRq7KskjCnKRbDCDiy7yG6GHHw59+/qNUd0IlTRQkItkU9Om8OqrMGKEj9DbtvV2f5EUKMhFsq24GG68Ef70J+8I7dAB7rxTHaFSaQpykVhOPNFvhJ5+Olx5JfToAR9+GLsqSSAFuUhMdev6VrgTJsDLL0PLlvD887GrkoRRkIvEtv2M0CVLoF49OPVUuOoq2Lw5dmWSEApykVxxzDEe5pdcAnfc4bsprlkTuypJAAW5SC6pWRPGj/fplrVrvb1/0iTQKYryHRTkIrmoVy+/EXrssfCTn8B558Gnn8auSnKUglwkVzVsCHPn+lLFadOgVStYtCh2VZKDFOQiuaxKFW8emj/fn59wAtx8szbfkm9RkIskQceOsGIF9Onjwd6tG2zYELsqyREKcpGk+N734NFH4aGHfHVLixbw9NOxq5IcoCAXSRIz6N/f92dp3Nhvig4bBl99FbkwiUlBLpJEhx8Or73mrf0TJkC7dvDWW7GrkkgU5CJJVb06jB4Ns2fDRx9BaSncc4/WnBcgBblI0nXvDitX+iZcQ4f6aUQbN8auSrJIQS6SDw44AGbO9Nb+mTN9861XXoldlWSJglwkXxQVwRVXeNPQPvtAly5w/fWwdWvsyiTDFOQi+aZNGz8btH9/7wrt1Mn3bZG8pSAXyUf77gsPPgiPPearWVq1gscfj12VZIiCXCSf9evnHaFHHum/HjjQj5eTvKIgF8l3TZrAvHne2v/QQ75MceXK2FVJGinIRQrB9gOf58717XDbtfNGIq05zwsKcpFCctJJvs/5SSd5a3/v3vDJJ7GrkhQpyEUKTUkJzJjha85nzPAboQsWxK5KUqAgFylE29ecv/YaVKsGnTvDTTdpn/OEUpCLFLLSUt9JsW9fGDkSTj4ZPvggdlVSQQpykUJXqxY88oivaFm82Nv7Z82KXZVUgIJcRP5/n/OlS6F+fejRw7fI3bw5dmWyF1IKcjO7wcxWmtkKM3vBzOqnqzARieCII3xUPmwY3HknHH88vPtu7KpkD1Idkd8eQmgRQmgFzAB+lXpJIhJVjRowbhw89RS8957v3fLYY7Grku+QUpCHED7b4ek+gLoLRPLFWWd5e3+LFnDeeTBggNr7c1TKc+RmdpOZrQPORyNykfxy8MHw8su+He6kSdC2rTcUSU7ZY5Cb2VwzW7WbR0+AEMKIEEIjYDJw6XdcZ7CZlZlZWXl5efr+BCKSWVWrwqhR8OKL8Nln0L49jB+v9v4cYiFN/zHM7BBgZgjhmD39bGlpaSgrK0vL54pIFpWX++qWWbOgVy944AHYf//YVRUMM1saQijd+fVUV60ctsPTM4H/SuV6IpLjSkrg2Wd9RcvMmd7e/+qrsasqeKnOkd+6bZplJXAK8LM01CQiuayoCH7xi2+39994o9r7I0p11UrvEMIx25Yg/jCEsCFdhYlIjtve3t+vn98M7dZN7f2RqLNTRCpvx/b+11/39v6ZM2NXVXAU5CKSmh3b+xs0gDPOUHt/linIRSQ9jjgCFi2CSy/1m6HHHaf2/ixRkItI+tSoAb/7HUyfDu+/D61bw+TJsavKewpyEUm/Xr28vb9VK7jgAvjJT+DzzyMXlb8U5CKSGQcfDC+9BL/6FTz8sK9yWbEidlV5SUEuIplTtSr85jfe3r9pE3To4Dsrqr0/rRTkIpJ5Xbr4aLxrV7jsMt9Z8ZNPYleVNxTkIpIdJSUwY4avaJk1y+fP58+PXVVeUJCLSPaYeXv/woXe3n/iiXDDDWrvT5GCXESyr21bb+8/91y/GdqtG2zQDh+VpSAXkThq1YI//MEPrHj9dZ9qee652FUlkoJcROIxg4su8vb++vXh9NNh+HDYsiV2ZYmiIBeR+La39w8ZArffDj/4AaxdG7uqxFCQi0huqFkT7r4bpk6F1at9quXJJ2NXlQgKchHJLX36wPLlcPjhcPbZMGwYfP117KpymoJcRHLPoYf6EXJXXgkTJnhH6Jo1savKWQpyEclN1arB6NHeRLR+PbRp46tcZBcKchHJbT16eHt/mzZw4YW+k+IXX8SuKqcoyEUk9zVsCH/6k58Nun0nxZUrY1eVMxTkIpIMVavCqFEwZw787/9C+/Zw773aSREFuYgkTdeuPtXSqZOvO+/XDz79NHZVUSnIRSR5DjjA2/lvucXXmrdpA0uWxK4qGgW5iCRTURFcey3Mmwdbt8Lxx8OYMQU51aIgF5FkO+44byA6/XS44go480zYuDF2VVmlIBeR5Nt/f5g+HcaOhRde8Pb+V1+NXVXWKMhFJD+Y+TFyCxdCjRp+aMVNNxXEoRUKchHJL23a+La455wDI0fCqafC3/4Wu6qMUpCLSP6pVQsmT4b774cFC6BlS19/nqcU5CKSn8xg4EBfllivHnTvDiNG+AqXPKMgF5H8dvTRHuYDB8LNN/vc+bp1satKKwW5iOS/f/s3uO8+ePRReOMNX9XyzDOxq0obBbmIFI5zz/U1540bQ8+e8POfwz/+EbuqlCnIRaSwNG0Kr70Gl18Od93lHaHvvRe7qpQoyEWk8FSv7iE+fTq8/z60bg1TpsSuqtIU5CJSuHr18p0Umzf3aZfBg+HLL2NXVWEKchEpbAcfDC+/7Btw3XcftGsHb78du6oKUZCLiBQX+5a4s2fDRx/5CUQPPpiYnRQV5CIi23Xv7ssTO3b0dec//jFs2hS7qj1KS5Cb2VVmFsysXjquJyISzUEH+Q6Ko0bBY49B27a+ZDGHpRzkZtYIOBn4S+rliIjkgCpV/KDnl16CL76ADh1g3LicnWpJx4h8DDAcyM0/oYhIZXXq5FMtJ5/sW+T27g1//3vsqnaRUpCb2ZnAhhDCG3vxs4PNrMzMysrLy1P5WBGR7KlXz9v577gDnn3W15wvXBi7qm/ZY5Cb2VwzW7WbR09gBPCrvfmgEMLEEEJpCKG0pKQk1bpFRLKnqMiPkVuwwH/9gx/AbbfBN9/ErgzYiyAPIXQLIRyz8wN4H2gCvGFma4GGwDIzOzCzJYuIRNKuHSxbBmedBddc4+eEfvRR7KoqP7USQngzhPD9EELjEEJjYD3QJoSQ30dxiEhhq10bpk6Fu+/2RqJWrfyfEWkduYhIRZnBkCGweLGfRtS1K/z619HOB01bkG8bmX+cruuJiOS8li2hrMwbh37zGw/0DRuyXoZG5CIiqdh3X5g0CR5+2EO9VSt47rmslqAgFxFJhwsv9CCvX99vgg4fDlu2ZOWjFeQiIulyxBGwaBEMHQq33+7LFNeuzfjHKshFRNKpZk2YMAGmTYPVq32q5cknM/qRCnIRkUw4+2w/tKJZM//1sGHw9dcZ+SgFuYhIpjRpAvPnw5VX+ii9QwdYsybtH6MgFxHJpGrVYPRomDEDysvh00/T/hFV035FERHZVY8e8N57UKNG2i+tEbmISLZkIMRBQS4ikngKchGRhFOQi4gknIJcRCThFOQiIgmnIBcRSTgFuYhIwlkIIfsfalYO/LmCb/sekO6WqHRcM5VrVOa9FXnP3v5sPUCHgmTm71i6ZLO2XP2upXqdir43E981SO37dkgIYdfT60MIiXgAE3PxmqlcozLvrch79vZngbLY/31z4ZGJv2NJrC1Xv2upXqei783Ed23bz6b9+5akqZVnc/SaqVyjMu+tyHsy8e8sn+Xyv69s1par37VUr1PR9ybmuxZlakVyi5mVhRBKY9chUggy8X1L0ohcMmdi7AJECkjav28akYuIJJxG5CIiCacgFxFJOAW5iEjCKchlF2Z2qJk9YGZPxK5FJJ+ZWS8zu8/MnjazUyp7HQV5gTCzB83sIzNbtdPrp5rZO2b2rpldCxBCeD+EMDBOpSLJVsHv2h9DCBcD/YG+lf1MBXnhmAScuuMLZlYFGA+cBhwFnGtmR2W/NJG8MomKf9dGbvv9SlGQF4gQwjzgk51ebge8u20EvhmYAvTMenEieaQi3zVzvwWeCyEsq+xnKsgLWwNg3Q7P1wMNzKyumd0DtDaz6+KUJpJXdvtdAy4DugFnm9mQyl68amq1ScLZbl4LIYSNQKX/UonILv7Vd20sMDbVi2tEXtjWA412eN4Q+CBSLSL5LKPfNQV5YVsCHGZmTcysGtAPeCZyTSL5KKPfNQV5gTCzx4CFQDMzW29mA0MIW4FLgeeB1cDUEMJbMesUSboY3zVtmiUiknAakYuIJJyCXEQk4RTkIiIJpyAXEUk4BbmISMIpyEVEEk5BLiKScApyEZGEU5CLiCTc/wFow6URxzRWlwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogx(TList, etaList, 'r-');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Here's a plot of our calculated result of $\\eta$ compared to the classical result.  As expected, the calculed $\\eta$ matches the classical result at the highest temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlrElEQVR4nO3deZxPZeP/8dc1ZuxLCred7DNjN4ylkiWp3JRSyDqWLBGSpSyFuw01qXAja0PE154HGnuWMdYxdkoh+07MjDm/P0716w4xM5+Z81nez8dj/piP+ZzzVj5vx3Wuc13GsixERMRz+TkdQEREUkZFLiLi4VTkIiIeTkUuIuLhVOQiIh5ORS4i4uH8nThprly5rKJFizpxahERj7Vt27ZzlmXl/vvrjhR50aJFiY6OduLUIiIeyxhz7G6va2hFRMTDqchFRDycilxExMOpyEVEPJyKXETEw6nIRUQ8nMcV+Y0bN5yOICLiVjyqyG/evElISAhdu3blypUrTscREXELHlXklmXxzDPPMGHCBIKDg1m6dKnTkUREHOdRRZ4pUyZGjx7Nxo0byZEjB40aNWLQoEFOxxIRcZQjj+inVGhoKNu3b+eDDz6gQYMGAMTHx+Pv748xxuF0IiJpyyOLHCB9+vQMHTr0z+/79OnDzz//zNixYylQoICDyURE0pZHDa38k2LFirFy5UqCgoKYOHEi2lRaRHyF1xR57969iYmJoUqVKnTu3Jl69erx448/Oh1LRCTVeU2RAxQvXpzIyEgmTpzI/v37SUxMdDqSiEiq86oiBzDG0LFjR3788UeKFy+OZVkMGDCAmJgYp6OJiKQKryvyP2TIkAGA48ePM3nyZCpXrszQoUO5deuWw8lERFzLa4v8D4UKFWLfvn00b96cYcOGUaVKFbZs2eJ0LBERl/H6Igd45JFHmDFjBkuXLuXy5cs0bdpUV+Yi4jV8osj/8OyzzxIbG8vChQvJkCEDCQkJ/PDDD07HEhFJEZ8qcoDs2bMTEhICwMSJE3nsscfo0KEDFy9edDiZiEjy+FyR/1W7du0YMGAA06ZNIygoiPnz5zsdSUQkyXy6yDNlysQHH3xAVFQUefPmpWnTpgwYMMDpWCIiSeKxa624UuXKlYmKimLUqFHUqVMHgLi4OAICArQIl4i4PZ++Iv+rgIAABg4cSPXq1QH7kf9nnnmGY8eOOZxMROSfqcjvISgoiA0bNlC2bFm+/PJLPe4vIm5LRX4P3bt3JzY2llq1avH666/zxBNPcPjwYadjiYjcwSVFboxpaIw5YIw5bIzxmruFRYoUYdmyZUybNo1jx47h56e/90TE/aS4mYwx6YAvgWeAIKCFMSYopcd1F8YY2rRpw9GjRylWrBiWZdG3b1+2b9/udDQREcA1V+TVgMOWZR21LCsO+AZo4oLjupWAgAAAfv31V2bOnEm1atUYMGAAv/32m8PJRMTXuaLICwC//OX747+/9j+MMZ2NMdHGmOizZ8+64LTOyJ8/P7GxsbRr146PPvqIihUrsn79eqdjiYgPc0WR322i9R37rFmWNcGyrBDLskJy587tgtM6J2fOnEyaNImVK1cSFxfHK6+8ws2bN52OJSI+yhVFfhwo9JfvCwInXXDcOy1ZAiNHQnQ03L6dKqdIivr167Nnzx6WLl1KxowZSUhIYPXq1U7HEhEf44oi3wqUNMY8aoxJDzQHFrnguHdavhz69YOqVeGRR6BJEwgPh127wKF53lmyZKFSpUoATJo0ibp169K6dWvOnTvnSB4R8T0pLnLLshKA14HlwD5gjmVZsSk97l19/jmcOAEREdCsGcTGQu/eULEi5MkDzZvDjBngUIm2b9+eIUOG8M033xAUFMTs2bOxrDtGmUREXMo4UTQhISFWdHS0aw72yy+wejWsWmVfsZ86BcZA9erQqBE89xyUL2+/lkZiYmLo0KEDW7dupU+fPowePTrNzi0i3ssYs82yrJA7Xvf4Iv+rxETYvh2WLrXH0/84R+HC8NJL9lV8aGialPrt27f57LPPqFWrFqGhody6dYv06dNrES4RSTbfKPK/O3UKvvsOFiywr9bj4qBQIbvQ07DUAbp168aBAweYMGECxYsXT5Nzioh3uVeRe/cz53nzQlgYLFoEZ87A9On2ePoXX0CNGlCkCPTpA1u3Qir/hVapUiWio6MpV64cn3zyCbfdYNaNiHgH774iv5fLl+1y//bb/3+lXqYMtG0LrVpBwYKpctrjx4/TtWtXlixZQrVq1Zg+fTqlS5dOlXOJiPfxzSvye8mRA1q3tsv89GmYOBFy5YKBA+3x9Keegq+/huvXXXraggULsmjRImbNmsWZM2dInz69S48vIr7JN6/I7+XIEXv64vTp8OOPkDWrPZb+2mtQrZpLx9MTEhLw9/fHsix69+5N8+bN/9zUQkTkbnRF/iCKF4d334XDh2HdOnjlFXv4pXp1qFQJxo2DK1dccip/f3uXvVOnTjFv3jxq1qxJnz59uO7ifwWIiPdTkd+Nnx88/jhMmgQnT8L48fZr3bpBvnzQseP/n9qYQvny5SM2NpYuXbrw6aefUq5cOSIjI11ybBHxDSry+8mWzR5a2bYNoqKgRQuYNcteJqBKFZgwIcVj6dmzZ2fs2LGsWbMGf39/WrVqpeVxReSBqcgflDF2ef9xlf7ll5CQYJd8wYLw1luQwo2aa9euza5du1i+fDmZMmUiPj6eFStWuOg3ICLeSkWeHDly2MMsO3fChg3QoAF8+ikUKwYvvghr1yZ7XnqmTJkoX748AJMnT+bpp5/m5Zdf5vTp0y78DYiIN1GRp4QxUKsWzJ5tz3Lp1w/WrIEnn7Rvjk6ZAilYpzwsLIz//Oc/LFy4kKCgIGbMmKFFuETkDipyVylUCD74AI4ft+el375tP1VauDCMGAEXLiT5kAEBAbz99tvs3LmT0qVL06ZNG3r16uX67CLi0VTkrpYpkz2rZfduiIy0x9UHD7YLvVevZI2jBwYGsn79esaMGUPLli0BuHnzJokOrcEuIu5FRZ5ajIG6de2VGHfvtsfOv/zSnqveqpW9GUYSpEuXjh49ehAaGgpA7969qV27NgcOHEiN9CLiQVTkaaFcOZg2DY4ehTfegIUL7cW7Gja011FPxrh39erViY2NpUKFCnz44YfEx8e7PreIeAQVeVoqVAhGj4aff4b337dnvdSrByEhMHdukrara9u2LXv37uXf//43AwcOJDQ0lH379qVedhFxWypyJ+TMaS/Q9dNP9o3Ra9fsNV3KlYOZMx94Y+m8efPy7bffMm/ePK5cuUKmTJlSN7eIuCUVuZMyZrRvjO7daz8tagy8+ioEBsLUqfCAwyVNmzblwIEDFC1aFMuy6NGjB+vXr0/d7CLiNlTk7iBdOnvj6N27Yd48e9XF9u2hVCl7CYBbtx7gEOkAOH36NEuWLOGJJ57g9ddf5+rVq6mdXkQcpiJ3J35+0LSpva7L4sWQJ4+9BECJEvauRg/wcFHevHmJiYnhjTfeYOzYsQQHB7Ns2bI0CC8iTlGRuyNjoFEj2LzZ3sGoaFHo0QMefRTCw+9b6FmzZiU8PJwffviBrFmzEhYWxo0bN9IkuoikPRW5OzPGXsdl3TpYvdoeO+/d256LPnbsfYdcatSowY4dO1i5ciWZM2cmPj6epUuX6jF/ES+jIvcExtjrt6xaZRd6sWLQvbs9hj5p0j/eFM2QIQNly5YFYOrUqTRq1IgXXniBkydPplF4EUltKnJP8+ST9hX68uWQNy906mRvHD19+n2nLbZv356RI0eyfPlygoKCmDRpkq7ORbyAitwT/THksnmzfVM0e3Zo2xaCg+Gbb+75YJG/vz99+/YlJiaGihUr0qlTJ7p3757G4UXE1VTknuyPm6LbttlPhvr72zsYVagA8+ff89H/EiVKsGrVKsaPH0+7du0AuHHjBrcf8EEkEXEvKnJv4OdnL8q1a5f9ZGhcnD2NsWZNe5OLu77Fj9dee41q1aoB0KdPH2rVqkVsbGxaJhcRF1CRe5N06ewr8thY+9H/X36xx9Sfeea+qy3Wrl2bI0eOUKlSJYYNG0ZcXFzaZBaRFFOReyN/f/vR/0OH4OOPYcsWe8eiVq3sFRjvokWLFuzdu5dmzZoxdOhQqlSpQkxMTBoHF5HkUJF7s0yZ7E2hjxyB/v3h//7PnuHSowfcZQ/Q3LlzExERweLFi0lISCBbtmwOhBaRpFKR+4KcOe1t6A4fttdwGTfOfqho6FC4cuWOH2/UqBGxsbF/LsLVrVs3Vq9e7UBwEXkQKSpyY0wzY0ysMSbRGBPiqlCSSvLnh//+1x5Df/ZZGDbMLvTw8DueEvXzs/9onD17lpUrV1K3bl06d+7MpUuX0j63iPyjlF6R7wGaAutckEXSSunSMGcOREXZUxV797aHXL755o4pi3ny5GHXrl289dZbfPXVVwQHB7No0SKHgovI3aSoyC3L2mdZljaN9FRVq8L338OKFZAjhz3jpXp12LDhf34sc+bMfPzxx2zZsoVcuXLRpUsXrl+/7lBoEfm7NBsjN8Z0NsZEG2Oiz549m1anlQfx1FP2Q0VTpsDx4/D44/a89EOH/ufHQkJCiI6OJjIykixZshAfH8/8+fP1mL+Iw+5b5MaY740xe+7y1SQpJ7Isa4JlWSGWZYXkzp07+YkldaRLB+3a2eU9fLi9lktQkL1Z9Llzf/5YQEAAgYGBAEyfPp2mTZvy3HPP8fPPPzsUXETuW+SWZdW3LKvsXb4WpkVASWOZM8OgQfYMlw4d7A0tSpSAkSPvWAe9Xbt2fPbZZ6xdu5bg4GDGjRtHYhI2kBYR19D0Q7m7vHlh/Hh7+7lataBfP3s99L/cEE2XLh09e/Zkz5491KhRg27dutG1a1eHg4v4npROP3zBGHMcqAEsNcYsd00scRvBwbB0Kaxc+b83RP+yufOjjz7K8uXLmTJlCh07dgTg+vXrJCQkOJVaxKekdNbKfMuyClqWlcGyrH9ZlvW0q4KJm6lf374hOnUqnDgBTzxhL8x15AgAxhjatWtH1apVAejduzehoaHs3LnTucwiPkJDK/Lg0qWz1z0/eNC+IbpihX1DtH//O54Qffrppzlx4gQhISG888473HyAjaNFJHlU5JJ0f9wQPXjQHmr5+GN727mvvvpzl6IXX3yRvXv30rp1a95//30qVqyoq3ORVKIil+TLn98eatm61Z7Z0rGj/ZDROvtB34cffpgpU6awfPlyAgICyJkzp7N5RbyUilxSLiTEvvk5a5Y957x2bWjWDH78EYAGDRqwe/duihQpgmVZdOrUieXLdV9cxFVU5OIaxkDz5rB/P7z3Hnz3nT1d8e234epVjDEAnDt3jg0bNtCwYUPatm3LhQsXHA4u4vlU5OJamTPDkCFw4IB9Vf7BB/b4+dSpkJhI7ty52bFjB4MGDWLmzJkEBgYyd+5cPeYvkgIqckkdBQvCjBmwaRMUKWKvgx4aCj/8QMaMGRk+fDjR0dEUKlSIN954Q4twiaSAilxSV/XqsHGjXeq//gqPPWbPdPnlFypUqMDmzZtZvXo1WbNmJS4ujjlz5ujqXCSJVOSS+vz87P1CDxywh10WLLDXRB8xAv+EBEqVKgVAREQEr7zyCg0aNODoPfYWFZE7qcgl7WTJYt8I3b/f3qFo8GB7CYBFi8CyaNu2LePGjWPLli2UK1eO8PBwbv8+L11E7k1FLmmvSBGYO9devyVDBmjSBJ57Dr/Dh+nSpQuxsbE8+eST9O7dm86dOzudVsTtqcjFOfXrw65d8Mkn9q5EZcvCgAEUypmTJUuWEBERQbdu3QC4du0acXFxDgcWcU8qcnFWQIC9Z+jBg9CyJXz0EZQujZk1i5YtWlClShXAXoQrJCSErVu3OhxYxP2oyMU95M1rzzXfuBHy5YNXX7WfEN21C4DGjRtz4cIFqlevzltvvcWNGzeczSviRlTk4l5q1IAtW2DCBNi7FypXhtdf59+1ahEbG0unTp0YNWoU5cuXZ9u2bU6nFXELKnJxP+nSQadO9v6h3brBuHFQqhQ5Zs9m/Jdfsnr1arJnz472fhWxqcjFfeXMCZ9/Dtu32+uev/YaVKvGkxkysG3bNgoXLoxlWYSFhbF48WKn04o4RkUu7q9CBVi7FmbOhFOnoGZNTLt2cOoU58+fJzo6msaNG9OyZUvOnj3rdFqRNKciF89gjP1o/4EDMGCAvWRu6dLkioggevNmhg0bxty5cwkMDCQiIkKP+YtPUZGLZ8ma1V5Rcc8e+8Zor16kr1mTwfXqsWPHDkqWLEm/fv20CJf4FBW5eKZSpWDZMvsJ0fPnoVYtgkePZsP8+axbt+7PRbgiIiJITEx0Oq1IqlKRi+cyBl58Efbtg379YMYM0gUFUXzlSrh9m1mzZtGqVSvq1KnDoUOHnE4rkmpU5OL5sma1nwjdtcu+Mdq1K1SvTpugIL766it27dpF+fLl+fjjj0lISHA6rYjLqcjFewQFwapVEBEBx49jQkMJ27qVvRs30rBhQ/r370+HDh2cTinicsaJu/shISFWdHR0mp9XfMjlyzB0qD0P/eGHsT76iHlZs1KsRAkqV67M1atXCQgIIGPGjE4nFXlgxphtlmWF/P11XZGLd8qRA8LD7YeJSpXCdOjAS2PGUDldOgD69OlDpUqV2Lhxo7M5RVxARS7erUIFWL8eJk+256BXqQK9evHSM89w48YNHnvsMXr27Mm1a9ecTiqSbCpy8X5+fvbmzwcO2Gu4jBnD0927s2fIELp368YXX3xB2bJliYqKcjqpSLKoyMV3PPywvQDXli1QsCDZOnbk8717WTdtGrlz5yZv3rxOJxRJFhW5+J6qVWHzZrvUd+7ksbAwourUofDDD2NZFm3atGHevHlOpxR5YCpy8U3p0kGXLvZwS+vWmJEjISiICzNmEBsby0svvcSLL77Ir7/+6nRSkftSkYtvy53bvhG6YQM89BCPtG3Llvz5+bB/f5YuXUpQUBBTpkzRIlzi1lI0j9wYMxL4NxAHHAHaW5Z16X7v0zxycUvx8fDZZ/b8c+Bg9+503LSJoz/+yL59+8iWLZvDAcXXpdY88pVAWcuyygMHgYEpPJ6IcwICoG9fe+2W+vUpNXIkay5eZEN4ONmyZSMuLo4pU6Zw+/Ztp5OK/I8UFbllWSssy/pj8YrNQMGURxJxWOHCsHAhLFiA35UrFG3WDDp1Ys7kyYSFhfH444+zd+9ep1OK/MmVY+RhwLJ7/aIxprMxJtoYE61dXMQjNGlibwDdty9MmcKrgwYxo3NnDh48SKVKlRgxYgTx8fFOpxS5/xi5MeZ74G4TbN+xLGvh7z/zDhACNLUeYNBdY+TicXbvtme5bNrEmZo16ZkjB7OXLePVV1/l66+/djqd+Ih7jZH73++NlmXVv8+B2wKNgHoPUuIiHql8eXtmy6RJ5Onfn2+uX6dFs2YU7t4dgCtXruDv70/mzJkdDiq+KEVDK8aYhkB/oLFlWTdcE0nETfn5QefO9tzz5s1p8u23VGrVCpYvp0+fPlSoUIE1a9Y4nVJ8UErHyL8AsgErjTE7jTHjXZBJxL3lyQPTp0NkpD3TpWFDXj1yBCshgTp16tClSxcuX77sdErxISmdtVLCsqxClmVV/P2ri6uCibi9unXtXYmGDaPOpk3sPn+eN+vWZeLEiQQHB7Np0yanE4qP0JOdIimRIQMMHgx79pC5enVGrVrFplKlKPLIIxQoUMDpdOIjVOQirlCiBCxfDt98Q7VLl/hhzx4KjxqFdekSLVu2ZNasWXrMX1KNilzEVYyBV16xnwzt2hW++IKLZcpwNDqali1b0rhxY44fP+50SvFCKnIRV3voIfjiC9iyhYfz5+eHQ4f4tEwZVkVGEhQUxH//+18SExOdTileREUuklqqVoWoKNJ99hm9Tpwg5vZtquXJw4gRI7h+/brT6cSLqMhFUpO/P/TsCfv2UaxRI1YeOcLGLFnIFhtLXFwcEydOJCEh4f7HEfkHKnKRtFCgAMybh1mwgELXr0PNmsx97jk6d+5MjRo12L17t9MJxYOpyEXS0h8LcfXsSYvISGY/9BDHDh6kSpUqDBkyhFu3bjmdUDyQilwkrWXLBuHhmKgoXi5alH1XrtAib16GDx9Ou3btnE4nHkhFLuKUkBDYupVHRo1i+oULLMuQgf4FCkBCApcvX+batWtOJxQPoSIXcZK/P7z5JsTG0rBePSqOHg2hofRp25Zy5cqxcuVKpxOKB1CRi7iDokVhyRKYPRtOnKD9okVkuHyZBg0aEBYWxsWLF51OKG5MRS7iLoyBl1+G/ft5rHNndl68yMDs2Zk+bRpBQUFs2LDB6YTiplTkIu7moYdg/HgybtjA+wULsjUxkdIJCRTJkMHpZOKmVOQi7qpWLdixg0ojRrDm6lUK1a+PNXYsLzdrxtSpU7UIl/xJRS7iztKnh3fegZgYCAnhcvfu/Lp8Oe3bt6dhw4b89NNPTicUN6AiF/EEJUvC99/z0PTprE2fni/9/Ni4Zg1ly5bl888/5/bt204nFAepyEU8hTHQujV++/fTrU0bYuPieDwxkZHDh3PjhrbM9WUqchFPkysXTJlC4VWr+K5AATafPUu2rl259csvjB07lvj4eKcTShpTkYt4qjp1MDEx5B88GObMYUFgIN27d6dq1aps27bN6XSShlTkIp4sY0YYNgx27uSVSpWYD5zZt4/Q0FAGDBjAb7/95nRCSQMqchFvEBQEa9fy/IQJ7M2YkfaWxUcffUTb1q2dTiZpQEUu4i38/KBTJx46cICJL7/M98Cg6GhYt45Lly5x5coVpxNKKlGRi3ibvHlh1izqLVtGeWOgdm36hIYSHBjI0qVLnU4nqUBFLuKtGjaEPXvgzTd57dAhcpw5Q6NGjWjVqhXnzp1zOp24kIpcxJtlyQKjRhG6dSvbg4N5F5gzcyaBpUuzZs0ah8OJq6jIRXxBlSqkj45m6Mcfsz19eipevkyJDRsgMdHpZOICKnIRX+HvD2+9RdnYWFbWqUPBwYOxatXixaeeYsKECSSq1D2WilzE1xQvDitWwLRpXNm/n4uRkbz22mvUq1OHw4cPO51OkkFFLuKLjIE2bchx8CCRLVowEdi+YQPlgoMZNWoUCQkJTieUJFCRi/iy3LkxERF0XLaMvfny0SAujjHvvcdvp045nUySQEUuItCwIQUOHGBB795svX6dbNWqcWv2bMLDw7l165bT6eQ+UlTkxpjhxpjdxpidxpgVxpj8rgomImksSxbMJ5/wr6go+Ne/WNy8Ob1796Zy+fJs3rzZ6XTyD1J6RT7SsqzylmVVBJYAQ1IeSUQcFRICUVG89OGHLA0I4OqhQ9SsWZM+vXtz/fp1p9PJXaSoyC3L+uviDVkAbSIo4g0CAqB/f57du5c9jz9OV8vi0/Bw2jz/vNPJ5C5SPEZujPmPMeYX4FV0RS7iXUqUIPuaNXw5ZQprs2bl3TVr4L33uHT6NBcvXnQ6nfzO3G8nbmPM90Deu/zSO5ZlLfzLzw0EMlqWNfQex+kMdAYoXLhwlWPHjiU7tIg44MwZ6NULZs0iLEcOlvn7M3biRF544QWnk/kMY8w2y7JC7nj9fkWehBMUAZZallX2fj8bEhJiRUdHu+S8IpLGvvuO7R060OHUKXYCLzVpwufjx5M3792u98SV7lXkKZ21UvIv3zYG9qfkeCLiAZ59lsqHDhHVowfvG8PihQsJKlmSyMhIp5P5rJSOkX9ojNljjNkNNADecEEmEXF3WbMSMGYMA7dsYWepUlS/do3So0bByZO46l/58uBSOmvlRcuyyv4+BfHflmWdcFUwEfEAVatSZs8evvvgAwquWYMVGMjzFSvyxeefaxGuNKQnO0UkZQICYMAA2L2bqxUqcGv3bnr07MkTVauyf79GW9OCilxEXKNkSbKvXcuyr75iWubM7Nu+nQply/L+8OHEx8c7nc6rqchFxHWMwYSF0eboUfY+/zzP377Nf4cP5+batU4n82oqchFxvX/9i3/Nn8/sJUvYljs32Ro04FbXrowcPpzffvvN6XReR0UuIqnnuefItX8/9OjBsvHj6TdkCBVLlGD9+vVOJ/MqKnIRSV3ZssFnn/H85s2sfPRR4k+e5IknnqB7u3ZcuXLl/u+X+1KRi0jaCA2l/oEDxAwdSq906Rg3bRpta9fWBtAuoCIXkbQTEECWd9/l07172VipEsN37oR69bgYHc25c+ecTuexVOQikvZKlaL6tm2UnTgRduzgzerVCSpalNkREXoyNBlU5CLiDGOgY0fYt49eTz5JkevXad6qFc/XqcOJE3pIPClU5CLirHz5KP/992yaM4dR2bKxcu1agooVY8WiRU4n8xgqchFxC/7NmvHmsWPsbtaMJ+PiCO7ZE1av1lDLA1CRi4j7yJmTEnPmsDAykgL+/iTWrUujIkUYPWIEt2/fdjqd21KRi4j7qVsXdu/m+htvkO6XX+g7eDA1ypRhz549TidzSypyEXFPmTOTLTychdHRzCpShB8PH6Zy+fIM7dOHW7duOZ3OrajIRcStmSpVaH7oEPsGD+ZlY5gWHk78pEmgsfM/qchFxP0FBJBr2DC+3r+f7TVqkPX117lVpw7vv/km169fdzqd41TkIuI5Spbk4fXrYcIElkdF8c4nn1CucGEiV6xwOpmjVOQi4ln8/KBTJxofPszaxx7D/8IF6j/9NB1feIFLly45nc4RKnIR8Uz58/PEunXsioigf+bMTF2wgNYhIeCD652ryEXEcxlDppYt+fD4caIaN+bDI0egQgUuLF7M6dOnnU6XZlTkIuL5cuak8sKFBH//Pdy+Td/GjQkqUoQZ48f7xJOhKnIR8R716kFMDG+1b0+ZW7do07Urz1apws8//+x0slSlIhcR75I5M4GTJ7M+Koox+fOzfscOgosXZ1lEhNPJUo2KXES8kl/VqvT46Sf2vPUWDRMTqdCtG0yejOWFOxKpyEXEewUEUPTjj/l23z7yV6xIYocOPJM7Nx/07Ut8fLzT6VxGRS4i3q9UKVi9mhvh4WS7coW3R4+mWtGi7Ni61elkLqEiFxHf4OdH1jfe4NuffmJe1ar8evIkVatVY2BYGDdv3nQ6XYqoyEXEtxQoQNMtW9g3eTJtMmbkmylTSBg0CDy4zFXkIuJ7jCFn+/ZMPnGCna++StbRo7lZrhzvhYVx9epVp9MlmYpcRHzXww+T4+uvYeVKIq9e5b0pUwjOn59lc+c6nSxJVOQiIvXr89yRI/zQogXZrl3j2WbNaFOnDufPn3c62QNxSZEbY/oaYyxjTC5XHE9EJM1lyUKNmTPZvmEDQ/LkYdaaNbQpWxY8YM2WFBe5MaYQ8BTg3c/AiohPyFCrFu8dP862Hj0Ydf48BAZy/vPPOXnihNPR7skVV+SfAv0A71+ZRkR8Q0AA5ceMIXD3bggO5q2ePQkqWpSvPvzQLRfhSlGRG2MaAycsy9rlojwiIu6jTBlYu5aB775LRcui48CB1C9ViqOHDjmd7H/ct8iNMd8bY/bc5asJ8A4w5EFOZIzpbIyJNsZEnz17NqW5RUTShp8fJYcOZdXRo4wvV46thw9TtkwZFo8Z43SyP5nk/jPBGFMOiARu/P5SQeAkUM2yrFP/9N6QkBArOjo6WecVEXGMZXF83Dj69enD6IQE8g0YQOLbb+OXOXOanN4Ys82yrJC/v57soRXLsmIsy8pjWVZRy7KKAseByvcrcRERj2UMBbt1Y+aJE+Rr1YrE//yHp3Pn5r2wMOLi4hyLpXnkIiJJ9cgjMHUqvy1cSB5jeHfKFKrky0fU6tWOxHFZkf9+ZX7OVccTEXF3WRo3JuL0aRY3acLFCxeoUbcufV94gRs3btz/zS6kK3IRkZTIkoVGCxYQGxlJp5w5+b8FC0hs3RrOnEmzCCpyEREXyFG3LuNPnWLXoEFkXbKEm2XKMKhxYy5dvJjq51aRi4i4Svr0ZBs+HHbuZE2+fHyweDHB+fKxaOLEVD2tilxExNUCA2kYE8OWfv3IFR9Pk86daV6pEmd+/TVVTqciFxFJDX5+hHz0EdGHDjGiVCnm79xJ28BA2LvX9ady+RFFRORPAcWK8c7+/ewcNYpP06eHW7dcfo5kP9mZEnqyU0R80q1bkCFDst/u8ic7RUQkiVJQ4v9ERS4i4uFU5CIiHk5FLiLi4VTkIiIeTkUuIuLhVOQiIh5ORS4i4uEceSDIGHMWOJbEt+UALrs4iiuOmZJjJOe9SXnPg/5sLkBryafOnzFXScts7vpZS+lxkvre1PisQco+b0Usy8p9x6uWZXnEFzDBHY+ZkmMk571Jec+D/iwQ7fT/X3f4So0/Y56YzV0/ayk9TlLfmxqftd9/1uWfN08aWlnspsdMyTGS896kvCc1/pt5M3f+75WW2dz1s5bS4yT1vR7zWXNkaEXcizEm2rrL+g0i4nqp8XnzpCtyST0TnA4g4kNc/nnTFbmIiIfTFbmIiIdTkYuIeDgVuYiIh1ORyx2MMcWMMV8ZY+Y6nUXEmxljnjfGTDTGLDTGNEjucVTkPsIYM9kYc8YYs+dvrzc0xhwwxhw2xgwAsCzrqGVZHZxJKuLZkvhZW2BZViegHfBKcs+pIvcdU4GGf33BGJMO+BJ4BggCWhhjgtI+mohXmUrSP2uDfv/1ZFGR+wjLstYBF/72cjXg8O9X4HHAN0CTNA8n4kWS8lkzto+AZZZlbU/uOVXkvq0A8Mtfvj8OFDDGPGKMGQ9UMsYMdCaaiFe562cN6AHUB14yxnRJ7sH9U5ZNPJy5y2uWZVnngWT/oRKRO9zrszYGGJPSg+uK3LcdBwr95fuCwEmHsoh4s1T9rKnIfdtWoKQx5lFjTHqgObDI4Uwi3ihVP2sqch9hjJkFbAJKG2OOG2M6WJaVALwOLAf2AXMsy4p1MqeIp3Pis6ZFs0REPJyuyEVEPJyKXETEw6nIRUQ8nIpcRMTDqchFRDycilxExMOpyEVEPJyKXETEw6nIRUQ83P8DasShx1gmAT4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Repeat the plot from above.\n",
    "plt.semilogx(TList, etaList, 'r-');\n",
    "\n",
    "# Add the classical result.\n",
    "nQ = (m*kB*np.array(TList)/(2*np.pi*hbar**2))**(3/2)\n",
    "etaClassical = np.log(n/nQ)\n",
    "plt.semilogx(TList, etaClassical, 'k--');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can also plot the temperature dependence of $\\mu$ form out calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeUklEQVR4nO3dfZzVY/7H8ddnZlZLK0lRpqEQbSLVdKNy052yqKybzbKSJTe5v6mt2Bu7blshqeRmY9308ysUS6VlWX6RiVaFMt1QhFg27Erp+v3xOa3Z2ampzpyuc873/Xw8zqM515mZ87mo9znne91ZCAEREUmWgtgFiIjI9qfwFxFJIIW/iEgCKfxFRBJI4S8ikkBFsQvYUvXr1w9NmjSJXYaISE6ZO3fuJyGEBpXbcyb8mzRpQllZWewyRERyipm9W1W7LvuIiCSQwl9EJIEU/iIiCaTwFxFJIIW/iEgCRQt/M+ttZovMrNzMfhGrDhGRJIoS/mZWCNwBHA20AE4xsxYxahERSaJY8/zbA+UhhKUAZjYJ6Au8WdNPdPvtsHYt7LUXHHAANGsGO+1U088iIpJbYoV/MbCiwv2VQIfK32Rmg4BBAHvttdc2PdGdd8LChf/ZVlICBx0EHTr4rV07qFdvm369iEhOihX+VkXbf50qE0KYAEwAKC0t3aZTZ+bPhzVrYPlyWLTou9vrr8PTT8PGs2z23x+6dYMePaBrV70YiEh+ixX+K4GSCvcbAx9k4onMYJddoFUrv1W0Zg2UlcGcOfDSS/DggzB+vP9M27b+QnDccdCxIxRoXpSI5BGLcYyjmRUBi4HuwPvAq8BPQwgLN/UzpaWlIdN7+6xbB6++CrNm+W32bFi/HvbYA/r2hX79/NNBrVoZLUNEpMaY2dwQQul/tcc6w9fMfgTcChQC94YQrt3c92+P8K/s88/90tDjj8NTT8GXX8IPfuCfBk49FY46Cr73ve1akojIVsm68N9aMcK/orVr4dln4bHHYMoU+PvfYbfd4OST4ac/hU6ddGlIRLLPpsJfcbWFatWCo4+GCRNg1Sp44gno2RMmToTDDoN99oGrroJly2JXKiJSPYX/NthhBzj2WHj4YfjoI/jjH6F5c7j+en8R6NkT/ud//NOCiEg2Uvinaeed4bTTYPp0n056zTXwzjvQvz8UF8Nll8GbNb50TUQkPQr/GlRSAldfDUuXwowZPjNozBg48EA48kgfK1i/PnaVIiIK/4woKPCZQI88Au+/DzfdBO++Cyee6JeFrr8ePvkkdpUikmQK/wxr0ACuvBLKy33K6P77w/Dh0LgxnHmmrzQWEdneFP7bSWGhLxSbNcv3GjrzTP9k0KYNHHEEPPkkbNgQu0oRSQqFfwQtWsDYsbByJYwa5QPFxx3nm81NnAjffBO7QhHJdwr/iOrWhUsv9UtCDzwARUUwcCA0bQojR8I//hG7QhHJVwr/LPC97/l2EfPm+SyhH/4QhgzxMwiGDIEPMrLlnYgkmcI/i5j5LKFZs3y30aOPhptv9hlCF1wAK1ZU/ztERLaEwj9LtW0LkybB4sXws5/5oTT77gvnnONjBCIi6VD4Z7l994W77vJxgbPO8gHhZs18tlB5eezqRCRXKfxzxN57+wyhpUvh/PN9X6EDDvBPBW+/Hbs6Eck1Cv8cU1wMt93mu4dedhk8+qhvH3H66f7CICKyJRT+OaphQ58Ounw5XH45TJ7snwTOPdfXD4iIbI7CP8c1aOB7By1Z4oPB994L++3n6wc+/jh2dSKSrRT+eaJRI99BdPFiP1ls9GifIjpiBHz2WezqRCTbKPzzTJMm/u7/zTf9wJnrrvMXgWuv9TOIRURA4Z+3DjjA1wnMm+fHTF51lU8RnTBBZwqIiMI/77VqBdOmwUsvfbdI7KCDYOpUCCF2dSISi8I/ITp1gr/+FR57zEO/Xz84/HB4+eXYlYlIDAr/BDHz0F+wAMaP97OGDz0UTjrJvxaR5FD4J1BRkV/+KS+HX/8ann7azxi44AJNDxVJCoV/gv3gB/CrX323b9D48T4ucMMN8PXXsasTkUxS+AsNG8K4cX45qGtXGDbMPwlMmaJBYZF8pfCXf2ve3GcGzZwJtWvDiSf6i4EOmRfJPwp/+S89e3rgjxvnh823bQtnnw0ffRS7MhGpKQp/qVJRkW8S9847vk/QxnMEbroJ1q6NXZ2IpEvhL5tVt64fJblwoV8CGjrUxwMefVTjASK5TOEvW2T//X1V8MyZsNNOcMIJ0KOH7yEkIrlH4S9bZeN4wNix/merVn6ewJo1sSsTka2h8JetVlQE553n20cPHAi33OIzhR58UJeCRHKFwl+2Wf36vkvoK69A48Zw2mlwxBHwxhuxKxOR6ij8JW3t2vkGcXfd5WMAbdrAxRfD55/HrkxENkXhLzWioMC3iFi82PcNGjPGzxSYOBE2bIhdnYhUpvCXGlWvHtxxB5SV+T5BAwdCly7w2muxKxORitIKfzM7ycwWmtkGMyut9NgwMys3s0Vm1qtCe1szm596bLSZWTo1SHZq3RpefNHf+S9Z4peGLrkEvvgidmUiAum/818A/Bh4oWKjmbUA+gMHAr2BsWZWmHp4HDAIaJa69U6zBslSBQUwYAAsWuSXgkaPhh/+UBvGiWSDtMI/hPBWCGFRFQ/1BSaFENaGEJYB5UB7M2sE1AkhzA4hBOB+oF86NUj2q1vX1wXMng0NGviGccceC8uXx65MJLkydc2/GFhR4f7KVFtx6uvK7ZIAHTrAq6/CqFHw/PO+TcSNN8K6dbErE0measPfzGaZ2YIqbn0392NVtIXNtG/quQeZWZmZla1evbq6UiUHFBX5RnFvvQW9e8MvfuHjAy+9FLsykWSpNvxDCD1CCC2ruE3dzI+tBEoq3G8MfJBqb1xF+6aee0IIoTSEUNqgQYPqSpUcUlLim8NNm+aDwF26+LbRn34auzKRZMjUZZ9pQH8zq2VmTfGB3TkhhFXAF2bWMTXL53Rgcy8ikueOO84Xhl15JfzhD75NxP33a0BYJNPSnep5vJmtBA4F/mRmMwBCCAuBR4A3genA4BDCt6kfOw+4Gx8EXgI8nU4Nkvtq1/ZzAl5/3c8MGDDAdwwtL49dmUj+spAjb7FKS0tDWVlZ7DIkwzZs8G0ihgyBb76Ba67xMYKiotiVieQmM5sbQiit3K4VvpJVCgp8TcDGAeEhQ6B9e60QFqlpCn/JSnvuCY895gvCVq3yF4ChQ+Gf/4xdmUh+UPhLVvvxj31AeOBAHxc4+GB49tnYVYnkPoW/ZL1dd/VxgGefBTPo3h1+/nP47LPYlYnkLoW/5IyuXf2gmKFD4b77fJ+gyZM1LVRkWyj8JafsuCPccINvE1FcDCedBMcfD++/H7sykdyi8Jec1Lq1Hx85ciTMnOn7BN1zjz4FiGwphb/krKIiuOIKmD/fXwzOOsunh773XuzKRLKfwl9y3r77+mDwHXf4BnEtW/rB8voUILJpCn/JCwUFcP75/imgXTtfKNazp84MENkUhb/klaZNYdYsGD/exwRatoRx43SIvEhlCn/JO2b+zn/BAujUyT8RdO8OS5fGrkwkeyj8JW/tvTfMmOELxObOhYMOgttv16cAEVD4S54z81lACxfC4YfDRRf5YjFtFy1Jp/CXRCgpgaee8gNj/vY33yPo1lv1KUCSS+EviWEGZ5zhnwK6dfNzArp1g2XLYlcmsv0p/CVxiovhiSd8RfBrr/mngLvu0roASRaFvySSGZx55nfrAgYNgmOPhQ8+iF2ZyPah8JdE23tvXxdw222+SrhlS5g0KXZVIpmn8JfEKyjwWUDz5sH++8Mpp8BPfgKffBK7MpHMUfiLpBxwALz4Ilx7rR8h2bIlPPlk7KpEMkPhL1JBUREMH+7nBey+Oxx3nJ8atmZN7MpEapbCX6QKrVr5C8CwYTBxos8Ieu652FWJ1ByFv8gm1KoF113nl4J22MHXBFxyCfzrX7ErE0mfwl+kGoce6oPBF17os4Jat4aysthViaRH4S+yBXbaCUaP9mmhX33lLwjXXQfffhu7MpFto/AX2Qrdu8Mbb8AJJ8CIEXDEEdoeQnKTwl9kK+26Kzz8MDzwgK8QbtUK7rtP20NIblH4i2wDMzj1VP8U0Lq1bxh38snw6aexKxPZMgp/kTTsvbdvC3HjjTB1qh8YM3Nm7KpEqqfwF0lTYSEMGeJnBtetC716aUqoZD+Fv0gNad3aj4vcOCW0XTs/OEYkGyn8RWrQjjv6lNDp0/36f7t2MHKkpoRK9lH4i2RAr14+E+jYY/2SUI8e8N57sasS+Y7CXyRD6teHKVPg3nt9RfDBB+usAMkeCn+RDDKDgQP92n+LFn5WwMCB8MUXsSuTpFP4i2wH++wDL7wAV18N998PbdpofyCJS+Evsp0UFcE11/jW0F9/7fsD3XQTbNgQuzJJorTC38xGmtnbZvaGmT1mZnUrPDbMzMrNbJGZ9arQ3tbM5qceG21mlk4NIrnm8MP9MlCfPjB0qA8Or1oVuypJmnTf+T8DtAwhHAwsBoYBmFkLoD9wINAbGGtmhamfGQcMApqlbr3TrEEk59SrB5Mnw4QJ8NJLPhj8xBOxq5IkSSv8QwgzQwjrU3dfBhqnvu4LTAohrA0hLAPKgfZm1gioE0KYHUIIwP1Av3RqEMlVZnD22fDaa9C4sX8SuPBCrQyW7aMmr/mfCTyd+roYWFHhsZWptuLU15Xbq2Rmg8yszMzKVq9eXYOlimSP5s3h5Zd9S4gxY6B9e1i4MHZVku+qDX8zm2VmC6q49a3wPSOA9cCDG5uq+FVhM+1VCiFMCCGUhhBKGzRoUF2pIjmrVi245RZ46in4+GMoLYVx47RNtGROUXXfEELosbnHzWwAcCzQPXUpB/wdfUmFb2sMfJBqb1xFu4gARx/t20QPGADnnw8zZsDdd/uCMZGalO5sn97AUKBPCOGfFR6aBvQ3s1pm1hQf2J0TQlgFfGFmHVOzfE4HpqZTg0i+2WMP/wQwapT/2aqVbxstUpPSveY/BtgZeMbM5pnZeIAQwkLgEeBNYDowOISwcWur84C78UHgJXw3TiAiKQUFcOmlvk30zjv73kDDh8O6dbErk3xhIUcuKpaWloYyLYmUBPrqK7j4YrjnHujUyY+Q3Guv2FVJrjCzuSGE0srtWuErkuVq1/br/g8/7DuFHnIIPP547Kok1yn8RXJE//6+JmCffeD44+Gii2Dt2thVSa5S+IvkkP328xXBl1wCt9/u+wO9807sqiQXKfxFcszGNQFTp8Ly5b5D6EMPxa5Kco3CXyRH9enjG8S1agWnngo//7kPDotsCYW/SA4rKYG//AVGjIA//MHPDF6wIHZVkgsU/iI5rqgIfvc7mDkT/v53fwG46y5tDSGbp/AXyRM9evhloMMOg0GD/MjIf/wjdlWSrRT+Inlkjz1g+nS47jo/L0DHRcqmKPxF8kxBAQwbBs8/79tBdOrks4N0GUgqUviL5KnOnWHePN8p9LLLoG9fHxMQAYW/SF6rV8+3grj1Vr8c1KYNzJkTuyrJBgp/kTxn5hvDvfii3+/SBUaP1mWgpFP4iyRE+/a+N1Dv3v5icNJJmg2UZAp/kQSpV8+3hRg50i8HtW0Lr78euyqJQeEvkjBmcMUVPhvo6699c7g779RloKRR+IskVOfO/q7/yCPh3HPhtNPgyy9jVyXbi8JfJMEaNPBzgn/3O5g0CUpL/cAYyX8Kf5GEKyjwjeFmzfIB4A4dYOLE2FVJpin8RQSArl39MlDHjjBwoN/++c/YVUmmKPxF5N8aNoRnnoGrr4b77vNPAW+/HbsqyQSFv4j8h8JCuOYaXxH84Yc+DqCTwvKPwl9EqnTUUb43UOvWflLYuef61FDJDwp/Edmk4mJ47jkYOtTXAnTp4ucGS+5T+IvIZhUVwQ03+Mrg8nLfHO6pp2JXJelS+IvIFunTB+bOhb32gmOOgV/+Er79NnZVsq0U/iKyxfbdF2bPhjPOgN/+1s8K+OST2FXJtlD4i8hW2XFHuPdePyT+hRf8MtArr8SuSraWwl9EtpoZnHUWvPSSTw097DAYO1abw+UShb+IbLO2bX0coGdPGDwYfvYz+Oqr2FXJllD4i0ha6tWDJ57wMYCHHvJVwYsWxa5KqqPwF5G0FRTAVVfBjBnw0UfQrh1MmRK7Ktkchb+I1JiePf2oyBYt4MQT/dCYdetiVyVVUfiLSI0qKfFZQBdcADffDN27w6pVsauSyhT+IlLjdtgBbr8dHnzQB4Rbt/ZjIyV7KPxFJGN++lOYMwd22cU/AYwcqemg2ULhLyIZdeCB8OqrcPzxMGQInHACrFkTuypR+ItIxtWpA4884mMA06bpkJhskFb4m9lvzewNM5tnZjPNbM8Kjw0zs3IzW2RmvSq0tzWz+anHRpuZpVODiOQGM7jsMj8p7NNPoX17eOyx2FUlV7rv/EeGEA4OIRwCPAn8EsDMWgD9gQOB3sBYMytM/cw4YBDQLHXrnWYNIpJDunb1QeDmzeHHP4bhw7U7aAxphX8IoeKVu9rAxqGcvsCkEMLaEMIyoBxob2aNgDohhNkhhADcD/RLpwYRyT0bp4OefTZcf73vDvrpp7GrSpa0r/mb2bVmtgI4ldQ7f6AYWFHh21am2opTX1du39TvHmRmZWZWtnr16nRLFZEs8v3vw4QJfnv+eT8r+PXXY1eVHNWGv5nNMrMFVdz6AoQQRoQQSoAHgQs2/lgVvypspr1KIYQJIYTSEEJpgwYNqu+NiOScs8+Gv/4V1q+HTp3gj3+MXVEyFFX3DSGEHlv4ux4C/gT8Cn9HX1LhscbAB6n2xlW0i0iCtW/v4wA/+QmcfrqvDbj5Zl8sJpmR7myfZhXu9gE2Tt6aBvQ3s1pm1hQf2J0TQlgFfGFmHVOzfE4HpqZTg4jkh91395lAl18OY8ZAt27aFiKT0r3mf0PqEtAbwFHAxQAhhIXAI8CbwHRgcAhh43j+ecDd+CDwEuDpNGsQkTxRVAS//z1MmuTX/9u08QNjpOZZyJG11qWlpaGsrCx2GSKynSxY4KuCly+HW2+F88/3tQKydcxsbgihtHK7VviKSFZq2dK3hejd23cIPeMM+Ne/YleVPxT+IpK16taFqVPhN7/xWUCdO/snAUmfwl9EslpBAfzyl35U5NKlfm7wM8/Erir3KfxFJCcccwyUlcGee/qloBtu0PbQ6VD4i0jO2G8/ePllOPlkGDYM+veHr76KXVVuUviLSE6pXRseeghuugkmT/ZVwcuWxa4q9yj8RSTnmMGVV8JTT8F77/m+QH/+c+yqcovCX0RyVq9ePh20USM46ii45RaNA2wphb+I5LT99oPZs6FfPz8sZsAArQfYEgp/Ecl5O+8M//u/8NvfwgMPwGGHwYoV1f9ckin8RSQvFBTAVVf5orDFi309wAsvxK4qeyn8RSSvHHecbwm9667QvTuMHatxgKoo/EUk7zRv7i8AvXrB4MF+YMzatbGryi4KfxHJS7vsAtOm+aWge+6BI4+ED3R01L8p/EUkbxUU+CDw5Mkwf76PA8yeHbuq7KDwF5G8d8IJvi3ETjvBEUfA3XfHrig+hb+IJMLG8wG6dvUxgMGD4ZtvYlcVj8JfRBKjXj3fEmLIEJ8F1KMHfPRR7KriUPiLSKIUFsKNN8LDD/sW0aWlMHdu7Kq2P4W/iCRS//7wf//ng8Jduvih8Umi8BeRxDrkEB8HaNcOTjkFRoyADRtiV7V9KPxFJNF23x1mzfJB4Ouu8w3i1qyJXVXmKfxFJPF22AHuvBPGjPEB4U6dYMmS2FVllsJfRAQ/IGbwYJg5E1atgvbt8/uAGIW/iEgF3br5vkANG/reQGPG5OfGcAp/EZFK9t3Xt4H40Y/gwgvhnHPyb0GYwl9EpAp16sDjj8Pw4XDXXb499Mcfx66q5ij8RUQ2oaAArr0WHnrIF4S1awfz5sWuqmYo/EVEqnHKKfDii/Dtt9C5M0yZErui9Cn8RUS2QNu2/u7/4IPhxBPh17/O7QVhCn8RkS3UsCH85S9wxhnwm9/ASSfBl1/GrmrbKPxFRLZCrVpw770wapQPCHfuDMuXx65q6yn8RUS2khlceqmvBn73XR8IfuGF2FVtHYW/iMg26tXLF4TttptPBZ0wIXZFW07hLyKShv339yMie/TwxWAXXADr1sWuqnoKfxGRNNWtC08+CZdfDnfc4SuDP/ssdlWbp/AXEakBhYXw+9/7YPDzz0PHjvDOO7Gr2rQaCX8zu8LMgpnVr9A2zMzKzWyRmfWq0N7WzOanHhttZlYTNYiIZIOBA3030E8/hQ4d4NlnY1dUtbTD38xKgJ7AexXaWgD9gQOB3sBYMytMPTwOGAQ0S916p1uDiEg2OewwHwhu1AiOOsrPCsg2NfHO/xZgCFBx09O+wKQQwtoQwjKgHGhvZo2AOiGE2SGEANwP9KuBGkREsso++/jOoL16wbnnwkUXwfr1sav6Tlrhb2Z9gPdDCH+r9FAxsKLC/ZWptuLU15XbN/X7B5lZmZmVrV69Op1SRUS2uzp1YNo0XxNw++1wzDHw+eexq3JF1X2Dmc0CGlbx0AhgOHBUVT9WRVvYTHuVQggTgAkApaWleXicgojku8JCXw3cogWcdx4ceig88QTst1/cuqp95x9C6BFCaFn5BiwFmgJ/M7PlQGPgNTNriL+jL6nwaxoDH6TaG1fRLiKS1846C555xs8E6NABnnsubj3bfNknhDA/hLB7CKFJCKEJHuxtQggfAtOA/mZWy8ya4gO7c0IIq4AvzKxjapbP6cDU9LshIpL9jjzSB4L32MMHgmOuCM7IPP8QwkLgEeBNYDowOITwberh84C78UHgJcDTmahBRCQbbTwisnt3XxF8ySVxBoIt5MjJxKWlpaGsrCx2GSIiNWL9erjiCrjtNujdGyZNgl12qfnnMbO5IYTSyu1a4SsiEkFREdx6K4wfD7Nm+UDwkiXb7/kV/iIiEZ1zDsycCR9+6APBzz+/fZ5X4S8iElnXrvDKK1C/vu8Oes89mX9Ohb+ISBZo1sy3hu7WzaeFXn65HxifKQp/EZEsUbcu/OlPfibAqFHQpw+sWZOZ51L4i4hkkaIi3wpi7FiYMQM6dYJVq2r+eRT+IiJZ6LzzYPp0Pylst91q/vdXu7ePiIjE0aOH3zJB7/xFRBJI4S8ikkAKfxGRBFL4i4gkkMJfRCSBFP4iIgmk8BcRSSCFv4hIAuXMYS5mthp4F6gPfBK5nFiS3HdIdv+T3HdIdv/T7fveIYQGlRtzJvw3MrOyqk6lSYIk9x2S3f8k9x2S3f9M9V2XfUREEkjhLyKSQLkY/hNiFxBRkvsOye5/kvsOye5/Rvqec9f8RUQkfbn4zl9ERNKk8BcRSaCcCX8z621mi8ys3Mx+EbueTDOzEjN7zszeMrOFZnZxqr2emT1jZu+k/tw1dq2ZYmaFZva6mT2Zup+kvtc1s8lm9nbq78ChSem/mV2a+ju/wMweNrPv53PfzexeM/vYzBZUaNtkf81sWCoHF5lZr2193pwIfzMrBO4AjgZaAKeYWYu4VWXceuDyEMIPgY7A4FSffwH8OYTQDPhz6n6+uhh4q8L9JPX9NmB6CKE50Ar/75D3/TezYuAioDSE0BIoBPqT332fCPSu1FZlf1MZ0B84MPUzY1P5uNVyIvyB9kB5CGFpCOEbYBLQN3JNGRVCWBVCeC319Rf4P/5ivN/3pb7tPqBflAIzzMwaA8cAd1doTkrf6wCHA/cAhBC+CSF8TkL6jx8vu6OZFQE7AR+Qx30PIbwA/L1S86b62xeYFEJYG0JYBpTj+bjVciX8i4EVFe6vTLUlgpk1AVoDrwB7hBBWgb9AALtHLC2TbgWGABsqtCWl7/sAq4E/pC573W1mtUlA/0MI7wO/B94DVgH/CCHMJAF9r2RT/a2xLMyV8Lcq2hIxR9XMfgBMAS4JIayJXc/2YGbHAh+HEObGriWSIqANMC6E0Br4ivy6zLFJqWvbfYGmwJ5AbTM7LW5VWaXGsjBXwn8lUFLhfmP8o2BeM7Pv4cH/YAjh0VTzR2bWKPV4I+DjWPVlUGegj5ktxy/xdTOzB0hG38H/vq8MIbySuj8ZfzFIQv97AMtCCKtDCOuAR4FOJKPvFW2qvzWWhbkS/q8CzcysqZntgA94TItcU0aZmeHXfN8KIYyq8NA0YEDq6wHA1O1dW6aFEIaFEBqHEJrg/6+fDSGcRgL6DhBC+BBYYWYHpJq6A2+SjP6/B3Q0s51S/wa64+NdSeh7RZvq7zSgv5nVMrOmQDNgzjY9QwghJ27Aj4DFwBJgROx6tkN/u+Af594A5qVuPwJ2w0f/30n9WS92rRn+73Ak8GTq68T0HTgEKEv9/38c2DUp/Qd+A7wNLAD+CNTK574DD+PjG+vwd/Y/31x/gRGpHFwEHL2tz6vtHUREEihXLvuIiEgNUviLiCSQwl9EJIEU/iIiCaTwFxFJIIW/iEgCKfxFRBLo/wEvH99GYp71ZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(TList, muList, 'b-');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}